/*
    Copyright 2007-2012 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Contains implementations and explicit instantiations of the
/// scatter_stage class. \ingroup mcrx

#include "config.h"
#include "mcrx-stage.h"
#include "mcrx-stage-impl.h"
#include "create_grid.h"
#include "arepo_grid.h"
#include "arepo_emission.h"

using namespace std;


template <template <typename> class grid_type> 
void mcrx::scatter_stage<grid_type>::setup_objects ()
{
  // Set up dust model (this one is painful)
  cout << "Setting up dust model" << endl;
  // First read wavelength vector
  CCfits::ExtHDU& lambda_hdu = open_HDU (*this->output_file_, "LAMBDA");
  read(lambda_hdu.column("lambda"),lambda_);

  // Read dust grains
  typename T_dust_model::T_scatterer_vector sv
    (read_dust_grains<typename T_dust_model::T_scatterer_vector>(this->p_, this->m_.units));
  this->model_.reset(new T_dust_model (sv));
  this->model_->set_wavelength(lambda_);

  // Now build grid
  cout << "Building grid" << endl;

  // The density generator is used to generate the densities of the
  // dust grains from the gas and metal densities in the simulation.  
  uniform_density_generator dg(this->p_);

  // by definition, we don't track intensities in this stage, so we
  // can set n_lambda to zero to save memory here
  const int decomp_level = this->p_.getValue("domain_level", int(1), true);
  const int n_lambda = 0; 
  T_grid* dummy=0;
  this->g_ = create_grid(*this->input_file_, dg, n_lambda, 
			 decomp_level, dummy).second;

  // Build emission
  cout << "Setting up emission" << endl;
  CCfits::ExtHDU& pdata_hdu = open_HDU(*this->input_file_, "PARTICLEDATA");
  this->emi_.reset(new T_emission(pdata_hdu, this->g_.get()));

// Build cameras based on -PARAMETERS HDUs in output file
  cout << "Setting up emergence" << endl;
  this->eme_.reset(new T_emergence (*this->output_file_));

  // if we are doing kinematic calculation, set the wavelength scale
  // of the cameras
  if(this->p_.getValue("use_kinematics", false, true)) {
    typename T_emergence::iterator stop= this->eme_->end();
    for (typename T_emergence::iterator c = this->eme_->begin();
	 c != stop; ++c)
      (*c)->set_dopplerscale(lambda_);
  }

  // tell cameras to allocate
  typename T_emergence::iterator stop= this->eme_->end();
  for (typename T_emergence::iterator c = this->eme_->begin();
       c != stop; ++c) {
    (*c)->allocate(this->emi_->zero_lambda());
    assert((*c)->get_image().size()>0);
  }
}


template <template <typename> class grid_type> 
void mcrx::scatter_stage<grid_type>::load_file ()
{
  // first set up general stuff
  setup_objects();

  // Only thing to load is camera images
  this->eme_->load_images(*this->output_file_, "");
}


template <template <typename> class grid_type> 
void mcrx::scatter_stage<grid_type>::load_dump (binifstream& dump_file)
{
  // first set up general stuff
  setup_objects();

  // Only thing to load is camera images
  this->eme_->load_dump(dump_file);
}


template <template <typename> class grid_type> 
void mcrx::scatter_stage<grid_type>::save_file ()
{
  // Get normalization and save images
  //const T_float normalization = 1./this->n_rays_;
  const T_float normalization = 1.;
  this->eme_->write_images(*this->output_file_, normalization, this->m_.units,
			   "",
			   this->p_.defined("compress_images")&&
			   this->p_.getValue("compress_images", bool ()),
			   this->p_.defined("write_images_parallel")&&
			   this->p_.getValue("write_images_parallel", bool ()));
}


template <template <typename> class grid_type> 
void mcrx::scatter_stage<grid_type>::save_dump (binofstream& dump_file) const
{
  this->eme_->write_dump(dump_file);
}


template <template <typename> class grid_type> 
bool mcrx::scatter_stage<grid_type>::shoot ()
{
  this->m_.mon_.set_images(*this->eme_);

  bool r= this->shoot_scatter (false);
  this->m_.mon_.clear_images();

  return r;
}

template <template <typename> class grid_type> 
void mcrx::scatter_stage<grid_type>::operator() ()
{
  // this is here so the function is instantiated
  this->run_stage();
}


// prevent instantiation of the xfer template in each stage. instead
// we do this in one place only
extern template class 
mcrx::xfer<mcrx::dust_model<mcrx::scatterer<mcrx::polychromatic_scatterer_policy,
					    mcrx::mcrx_rng_policy>,
			    mcrx::cumulative_sampling,
			    mcrx::mcrx_rng_policy>,
	   mcrx::dummy_grid>;
extern template class 
mcrx::xfer<mcrx::dust_model<mcrx::scatterer<mcrx::polychromatic_scatterer_policy,
					    mcrx::mcrx_rng_policy>,
			    mcrx::cumulative_sampling,
			    mcrx::mcrx_rng_policy>,
	   mcrx::T_full_sed_adaptive_grid>;
#ifdef WITH_AREPO
extern template class 
mcrx::xfer<mcrx::dust_model<mcrx::scatterer<mcrx::polychromatic_scatterer_policy,
					    mcrx::mcrx_rng_policy>,
			    mcrx::cumulative_sampling,
			    mcrx::mcrx_rng_policy>,
	   mcrx::T_full_sed_arepo_grid>;
#endif

// explicit instantiations
template class mcrx::scatter_stage<mcrx::adaptive_grid>;

#ifdef WITH_AREPO
template class mcrx::scatter_stage<mcrx::arepo_grid>;
#endif
